/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
**                                                                    **
** (C) Copyright 1999, The Regents of the University of California    **
** All Rights Reserved.                                               **
**                                                                    **
** Commercial use of this program without express permission of the   **
** University of California, Berkeley, is strictly prohibited.  See   **
** file 'COPYRIGHT'  in main directory for information on usage and   **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
**                                                                    **
** Developed by:                                                      **
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
**                                                                    **
** ****************************************************************** */
                                                                        
// Original implementation: José Abell (UANDES), Massimo Petracca (ASDEA)
//
// ASDPlasticMaterial3D
//
// Fully general templated material class for plasticity modeling

#ifndef RoundedMohrCoulomb_YF_H
#define RoundedMohrCoulomb_YF_H

#include "../YieldFunctionBase.h"
#include "cmath"
#include <iostream>
#include "../AllASDModelParameterTypes.h"


template<class EtaHardeningType>
class RoundedMohrCoulomb_YF : public YieldFunctionBase<RoundedMohrCoulomb_YF<EtaHardeningType>> // CRTP
{
public:

    static constexpr const char* NAME = "RoundedMohrCoulomb_YF";


    RoundedMohrCoulomb_YF( ):
        YieldFunctionBase<RoundedMohrCoulomb_YF<EtaHardeningType>>::YieldFunctionBase() 
        {}

    YIELD_FUNCTION 
    {


        double p = -sigma.meanStress();
        double q = sigma.stressDeviatorQ();
        double theta = sigma.lodeAngle();

        double eta = GET_TRIAL_INTERNAL_VARIABLE(EtaHardeningType).value();

        double m = GET_PARAMETER_VALUE(RMC_m);
        double qa = GET_PARAMETER_VALUE(RMC_qa);
        double pc = GET_PARAMETER_VALUE(RMC_pc);
        double e = GET_PARAMETER_VALUE(RMC_e);


        double s0 = sigma(0);
        double s1 = sigma(1);
        double s2 = sigma(2);
        double s3 = sigma(3);
        double s4 = sigma(4);
        double s5 = sigma(5);

        // double yf = q * pow(1 + q / qa, m)  - eta * (p - pc) / g(theta, e);
        double F_voigt = q * pow(1 + q / qa, m) * g(theta, e)  - eta * (p + pc) ;
        // double F_voigt = -eta*(-pc - 1.0/3.0*s0 - 1.0/3.0*s1 - 1.0/3.0*s2) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))));

        if(F_voigt != F_voigt)
        {
            std::cout << "RoundedMohrCoulomb_YF - NaN detected! \n";
            std::cout << "  sigma = " << sigma.transpose() << " \n";
            std::cout << "  eta = " << eta << " \n";
            std::cout << "  m = " << m << " \n";
            std::cout << "  qa = " << qa << " \n";
            std::cout << "  pc = " << pc << " \n";
            std::cout << "  e = " << e << " \n";
            std::cout << "  F_voigt = " << F_voigt << " \n";
            std::cout << "  p = " << p << " \n";
            std::cout << "  q = " << q << " \n";
            std::cout << "  theta = " << theta << " \n";
            std::cout << "  g(theta, e) = " << g(theta, e) << " \n";
        }

        return F_voigt;
    }

    YIELD_FUNCTION_STRESS_DERIVATIVE
    {  

        double eta = GET_TRIAL_INTERNAL_VARIABLE(EtaHardeningType).value();

        double m = GET_PARAMETER_VALUE(RMC_m);
        double qa = GET_PARAMETER_VALUE(RMC_qa);
        double pc = GET_PARAMETER_VALUE(RMC_pc);
        double e = GET_PARAMETER_VALUE(RMC_e);
        
        double s0 = sigma(0);
        double s1 = sigma(1);
        double s2 = sigma(2);
        double s3 = sigma(3);
        double s4 = sigma(4);
        double s5 = sigma(5);

        vv_out(0) = (1.0/3.0)*eta + (1.0/2.0)*m*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(2*s0 - s1 - s2)/(qa*(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa)*((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(2*s0 - s1 - s2)/(((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(-1.0/36.0*sqrt(6)*(2 - 2*pow(e, 2))*(-1.0/432.0*(12*s0 - 6*s1 - 6*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - 1.0/2.0*(12*s1*s2 - 12*pow(s3, 2) - 2*pow(s0 - s1, 2) - 2*pow(s0 - s2, 2) - 2*pow(s1 - s2, 2) + (4.0/3.0)*pow(s0 + s1 + s2, 2) - 2*(s0 + s1 + s2)*(4*s0 - 2*s1 - 2*s2))*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 3.0/2.0) - 1.0/72.0*sqrt(6)*(2 - 2*pow(e, 2))*(12*s0 - 6*s1 - 6*s2)*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/(sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))*(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) - (2*e - 1)*((1.0/432.0)*(4 - 4*pow(e, 2))*(-1.0/216.0*(12*s0 - 6*s1 - 6*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (12*s1*s2 - 12*pow(s3, 2) - 2*pow(s0 - s1, 2) - 2*pow(s0 - s2, 2) - 2*pow(s1 - s2, 2) + (4.0/3.0)*pow(s0 + s1 + s2, 2) - 2*(s0 + s1 + s2)*(4*s0 - 2*s1 - 2*s2))*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2) + (1.0/432.0)*(4 - 4*pow(e, 2))*(12*s0 - 6*s1 - 6*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)))/sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/pow((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))), 2) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*(-1.0/216.0*(12*s0 - 6*s1 - 6*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (12*s1*s2 - 12*pow(s3, 2) - 2*pow(s0 - s1, 2) - 2*pow(s0 - s2, 2) - 2*pow(s1 - s2, 2) + (4.0/3.0)*pow(s0 + s1 + s2, 2) - 2*(s0 + s1 + s2)*(4*s0 - 2*s1 - 2*s2))*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2) + (1.0/216.0)*(4 - 4*pow(e, 2))*(12*s0 - 6*s1 - 6*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))));
        vv_out(1) = (1.0/3.0)*eta + (1.0/2.0)*m*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(-s0 + 2*s1 - s2)/(qa*(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa)*((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(-s0 + 2*s1 - s2)/(((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(-1.0/36.0*sqrt(6)*(2 - 2*pow(e, 2))*(-1.0/432.0*(-6*s0 + 12*s1 - 6*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - 1.0/2.0*(12*s0*s2 - 12*pow(s4, 2) - 2*pow(s0 - s1, 2) - 2*pow(s0 - s2, 2) - 2*pow(s1 - s2, 2) - 2*(-2*s0 + 4*s1 - 2*s2)*(s0 + s1 + s2) + (4.0/3.0)*pow(s0 + s1 + s2, 2))*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 3.0/2.0) - 1.0/72.0*sqrt(6)*(2 - 2*pow(e, 2))*(-6*s0 + 12*s1 - 6*s2)*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/(sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))*(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) - (2*e - 1)*((1.0/432.0)*(4 - 4*pow(e, 2))*(-1.0/216.0*(-6*s0 + 12*s1 - 6*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (12*s0*s2 - 12*pow(s4, 2) - 2*pow(s0 - s1, 2) - 2*pow(s0 - s2, 2) - 2*pow(s1 - s2, 2) - 2*(-2*s0 + 4*s1 - 2*s2)*(s0 + s1 + s2) + (4.0/3.0)*pow(s0 + s1 + s2, 2))*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2) + (1.0/432.0)*(4 - 4*pow(e, 2))*(-6*s0 + 12*s1 - 6*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)))/sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/pow((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))), 2) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*(-1.0/216.0*(-6*s0 + 12*s1 - 6*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (12*s0*s2 - 12*pow(s4, 2) - 2*pow(s0 - s1, 2) - 2*pow(s0 - s2, 2) - 2*pow(s1 - s2, 2) - 2*(-2*s0 + 4*s1 - 2*s2)*(s0 + s1 + s2) + (4.0/3.0)*pow(s0 + s1 + s2, 2))*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2) + (1.0/216.0)*(4 - 4*pow(e, 2))*(-6*s0 + 12*s1 - 6*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))));
        vv_out(2) = (1.0/3.0)*eta + (1.0/2.0)*m*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(-s0 - s1 + 2*s2)/(qa*(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa)*((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(-s0 - s1 + 2*s2)/(((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(-1.0/36.0*sqrt(6)*(2 - 2*pow(e, 2))*(-1.0/432.0*(-6*s0 - 6*s1 + 12*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - 1.0/2.0*(12*s0*s1 - 12*pow(s5, 2) - 2*pow(s0 - s1, 2) - 2*pow(s0 - s2, 2) - 2*pow(s1 - s2, 2) - 2*(-2*s0 - 2*s1 + 4*s2)*(s0 + s1 + s2) + (4.0/3.0)*pow(s0 + s1 + s2, 2))*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 3.0/2.0) - 1.0/72.0*sqrt(6)*(2 - 2*pow(e, 2))*(-6*s0 - 6*s1 + 12*s2)*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/(sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))*(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) - (2*e - 1)*((1.0/432.0)*(4 - 4*pow(e, 2))*(-1.0/216.0*(-6*s0 - 6*s1 + 12*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (12*s0*s1 - 12*pow(s5, 2) - 2*pow(s0 - s1, 2) - 2*pow(s0 - s2, 2) - 2*pow(s1 - s2, 2) - 2*(-2*s0 - 2*s1 + 4*s2)*(s0 + s1 + s2) + (4.0/3.0)*pow(s0 + s1 + s2, 2))*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2) + (1.0/432.0)*(4 - 4*pow(e, 2))*(-6*s0 - 6*s1 + 12*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)))/sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/pow((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))), 2) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*(-1.0/216.0*(-6*s0 - 6*s1 + 12*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (12*s0*s1 - 12*pow(s5, 2) - 2*pow(s0 - s1, 2) - 2*pow(s0 - s2, 2) - 2*pow(s1 - s2, 2) - 2*(-2*s0 - 2*s1 + 4*s2)*(s0 + s1 + s2) + (4.0/3.0)*pow(s0 + s1 + s2, 2))*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2) + (1.0/216.0)*(4 - 4*pow(e, 2))*(-6*s0 - 6*s1 + 12*s2)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))));
        vv_out(3) = 3*m*s3*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))/(qa*(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa)*((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))) + 3*M_SQRT2*s3*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))/(((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(-1.0/2.0*sqrt(6)*s3*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/(sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))*(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) - 1.0/36.0*sqrt(6)*(2 - 2*pow(e, 2))*(-1.0/12.0*s3*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - 1.0/2.0*(-24*s0*s3 + 24*s4*s5)*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 3.0/2.0) - (2*e - 1)*((1.0/12.0)*s3*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (1.0/432.0)*(4 - 4*pow(e, 2))*(-1.0/6.0*s3*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (-24*s0*s3 + 24*s4*s5)*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2))/sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/pow((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))), 2) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/6.0)*s3*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (1.0/216.0)*(4 - 4*pow(e, 2))*(-1.0/6.0*s3*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (-24*s0*s3 + 24*s4*s5)*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))));
        vv_out(4) = 3*m*s4*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))/(qa*(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa)*((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))) + 3*M_SQRT2*s4*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))/(((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(-1.0/2.0*sqrt(6)*s4*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/(sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))*(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) - 1.0/36.0*sqrt(6)*(2 - 2*pow(e, 2))*(-1.0/12.0*s4*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - 1.0/2.0*(-24*s1*s4 + 24*s3*s5)*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 3.0/2.0) - (2*e - 1)*((1.0/12.0)*s4*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (1.0/432.0)*(4 - 4*pow(e, 2))*(-1.0/6.0*s4*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (-24*s1*s4 + 24*s3*s5)*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2))/sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/pow((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))), 2) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/6.0)*s4*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (1.0/216.0)*(4 - 4*pow(e, 2))*(-1.0/6.0*s4*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (-24*s1*s4 + 24*s3*s5)*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))));
        vv_out(5) = 3*m*s5*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))/(qa*(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa)*((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))) + 3*M_SQRT2*s5*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))/(((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + pow(2*e - 1, 2))*(-1.0/2.0*sqrt(6)*s5*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/(sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))*(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))) - 1.0/36.0*sqrt(6)*(2 - 2*pow(e, 2))*(-1.0/12.0*s5*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - 1.0/2.0*(-24*s2*s5 + 24*s3*s4)*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 3.0/2.0) - (2*e - 1)*((1.0/12.0)*s5*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (1.0/432.0)*(4 - 4*pow(e, 2))*(-1.0/6.0*s5*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (-24*s2*s5 + 24*s3*s4)*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2))/sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/pow((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))), 2) + (1.0/2.0)*M_SQRT2*pow(1 + (1.0/2.0)*M_SQRT2*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/qa, m)*((1.0/6.0)*s5*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (1.0/216.0)*(4 - 4*pow(e, 2))*(-1.0/6.0*s5*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 2) - (-24*s2*s5 + 24*s3*s4)*(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/pow((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2), 2))*sqrt(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2))/((1.0/36.0)*sqrt(6)*(2 - 2*pow(e, 2))*sqrt(pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3))/sqrt((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2)) + (2*e - 1)*sqrt(5*pow(e, 2) - 4*e + (1.0/216.0)*(4 - 4*pow(e, 2))*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3)/((1.0/216.0)*pow(6*pow(s3, 2) + 6*pow(s4, 2) + 6*pow(s5, 2) + pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2), 3) + pow(6*s0*s1*s2 - 6*s0*pow(s3, 2) - 6*s1*pow(s4, 2) - 6*s2*pow(s5, 2) + 12*s3*s4*s5 + (2.0/9.0)*pow(s0 + s1 + s2, 3) - (s0 + s1 + s2)*(pow(s0 - s1, 2) + pow(s0 - s2, 2) + pow(s1 - s2, 2)), 2))));



        return vv_out;
    }

    YIELD_FUNCTION_XI_STAR_H_STAR
    {
        double dbl_result = 0.0;
        double eta = GET_TRIAL_INTERNAL_VARIABLE(EtaHardeningType).value();
        double pc = GET_PARAMETER_VALUE(RMC_pc);
        double s0 = sigma(0);
        double s1 = sigma(1);
        double s2 = sigma(2);
        double dFd_eta = pc + (1.0/3.0)*s0 + (1.0/3.0)*s1 + (1.0/3.0)*s2;
        dbl_result = dFd_eta*GET_INTERNAL_VARIABLE_HARDENING(EtaHardeningType).value();
        return dbl_result;
    }

  
    using internal_variables_t = std::tuple<EtaHardeningType>;

    using parameters_t = std::tuple<RMC_m,RMC_qa,RMC_pc,RMC_e>;

private:

    double g(double lode, double e) const // Willam-Warnke function
    {
        double coslode = cos(lode);
        double coslode2 = coslode * coslode;
        double num = 4 * (1 - e * e) * coslode2  + (2 * e - 1) * (2 * e - 1);
        double delta = 4 * (1. - e * e) * coslode2 + 5 * e * e - 4 * e;
        if (delta < 0)
        {
            delta = 0;
        }
        double den = 2 * (1 - e * e) * coslode + (2 * e - 1) * sqrt(delta);
        if(abs(den)<1e-20)
        {
            return 1.;
        }
        else
        {
            return num / den;
        }
    }

    double dg_dlode(double theta, double e) const //Derivative of above.
    {
        double dg = (2 * pow(e, 2) - 2) * (4 * ((2 * e - 1) * sqrt(5 * pow(e, 2) - 4 * e + (-4 * pow(e, 2) + 4) * pow(cos(theta), 2)) + (-2 * pow(e, 2) + 2) * cos(theta)) * sqrt(5 * pow(e, 2) - 4 * e + (-4 * pow(e, 2) + 4) * pow(cos(theta), 2)) * cos(theta) - ((4 * e - 2) * cos(theta) + sqrt(5 * pow(e, 2) - 4 * e + (-4 * pow(e, 2) + 4) * pow(cos(theta), 2))) * (pow(2 * e - 1, 2) + (-4 * pow(e, 2) + 4) * pow(cos(theta), 2))) * sin(theta) / (pow((2 * e - 1) * sqrt(5 * pow(e, 2) - 4 * e + (-4 * pow(e, 2) + 4) * pow(cos(theta), 2)) + (-2 * pow(e, 2) + 2) * cos(theta), 2) * sqrt(5 * pow(e, 2) - 4 * e + (-4 * pow(e, 2) + 4) * pow(cos(theta), 2)));
        return dg;
    }

    static VoigtVector vv_out; //For returning VoigtVector's
};

template <class EtaHardeningType>
VoigtVector RoundedMohrCoulomb_YF<EtaHardeningType>::vv_out;

//Declares this YF as featuring an apex
template<class EtaHardeningType>
struct yf_has_apex<RoundedMohrCoulomb_YF<EtaHardeningType>> : std::true_type {};

#endif



